Model Hamiltonian for Topological Insulators 
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In this paper we give the full microscopic derivation of the model Hamiltonian for the three 
dimensional topological insulators in the Bi2Se3 family of materials {Bi2Se3, Bi2Te3 and Sb2Te3). 
We first give a physical picture to understand the electronic structure by analyzing atomic orbitals 
and applying symmetry principles. Subsequently, we give the full microscopic derivation of the 
model Hamiltonian introduced by Zhang et al [1] based both on symmetry principles and the k • 
p perturbation theory. Two different types of k^ terms, which break the in-plane full rotation 
symmetry down to three fold rotation symmetry, are taken into account. Effective Hamiltonian is 
derived for the topological surface states. Both the bulk and the surface models are investigated in 
the presence of an external magnetic field, and the associated Landau level structure is presented. For 
more quantitative fitting to the first principle calculations, we also present a new model Hamiltonian 
including eight energy bands. 

PACS numbers: 71.15.-m, 71.18.+y, 73.20.-r, 73.61. Le 



I. INTRODUCTION 



Recently, topological insulators (TI) have been in- 
vestigated intensively both theoretically and experimen- 
tally. J2|-[3 These insulators are fully gapped in the 
bulk, but have gapless edge or surface states which are 
topologically protected by the time reversal symmetry. 
Topological insulator was first theoretically predicted [5| 
and experimentally observed [6] in the HgTe quantum 
wells. Transport measurements [6, 7] show the existence 
of the gapless edge channel, which demonstrates that 
HgTe/CdTe quantum well is a two-dimensional (2D) TI 
with quantum spin Hall effect. Later, Bi^Sbi-x was sug- 
gested to be a three-dimensional (3D) TI[8] with topolog- 
ically non- trivial surface states, which were observed W 
angle-resolved photoemission spectroscopy (ARPES)[9|. 
However, Bi^Sbi-x has a small energy gap, alloy disor- 
der and rather complicated surface states. More recently, 
new TIs with large bulk gaps ~ O.SeV and single Dirac 
cone surface states have been theoretically predicted 
for 5^2^63 [IJ^^sTes PL] and ^^2^63 P, M. ARPES 
measurement flOl UA\ indeed shows the single Dirac cone 
with linear dispersion around the F point in both Bi2Ses 
and Bi2Tes. Current research on these materials is de- 
veloping rapidly. |12l-l3Q| 



For deeper understanding and quantitative predictions 
of novel phenomena associated with the TIs, it is highly 
desirable to construct standard models for both 2D and 
3D TIs. Bernevig, Hughes and Zhang (BHZ)[5] con- 
structed the model Hamiltonian for the 2D TI in HgTe 
quantum wells. This model Hamiltonian demonstrates 
the basic mechanism of TI behavior through band inver- 
sion induced by spin-orbit coupling (SOC). It has been 
applied successfully for quantitative predictions of the he- 
lical edge states and properties under magnetic fields. |3l| 



Zhang et al[l] derived a model Hamiltonian for the 3D 
TI Bi2Ses, Bi2Tes and Sb2Tes, and obtained topolog- 
ical surface states consisting of a single Dirac cone. In- 
terestingly, in the thin film limit, the 3D TI model re- 
duces exactly to the 2D TI model by BHZ[32-'3^. In 
this paper, we give the full microscopic derivation of 
our model Hamiltonian, first by constraining its form by 
symmetry principles and a careful analysis of the rele- 
vant atomic orbitals. Subsequently, we determine the 
parameters of our model Hamiltonian by a systematic 
k • p expansion near the F point, and comparison with 
the ab initio calculations [1]. Furthermore the higher or- 
der k^ terms neglected in Ref. [l|, are also included in 
the derivation in order to recover the crystal C3 rotation 
symmetry |35|. Compared to the symmetry arguments 
given in Ref.fl], the new derivation given in this paper 
determines all the parameters of our model Hamiltonian 
by the wavefunctions from ab initio calculation, so that 
no fitting is required and no ambiguity is introduced. 
As an application of our model Hamiltonian, we study 
the bulk and surface Landau level spectra in a magnetic 
field. The surface Landau levels have V^ field depen- 
dence, as is expected from the Dirac-type dispersion of 
the surface states. The gap of 0th Landau level can be 
as large as bOmeV for lOT magnetic field, which suggests 
that the topological magneto-electric effect [36|, [S^ can be 
observable at such energy scales. Furthermore, we pro- 
pose a more quantitative description of the Bi2Ses fam- 
ily of TIs by going beyond the four bands and present 
a new model Hamiltonian with eight bands. Recently, 
our model Hamiltonian has been applied successfully for 
understand a number of experiments, including the STM 
study of the topological surface states [14, 20], STM study 
of the surface bound states [30], STM study of the quasi- 
particle interference [14, 20, 38], crossover from 3D to 2D 
topological insulators [21.. ■32.-.34] . and the Landau level of 



the topological surface states [22|, |23|. 

The paper is organized as follows. In sec. [Ill we first 
present the lattice structure and the symmetry proper- 
ties of Bi2Se3 crystal. Then we turn to the electronic 
band structure and discuss about the atomic orbital pic- 
ture, which is helpful to capture the essential physics in 
the long wave length limit. Keeping such atomic orbital 
picture in mind, we investigate in detail the properties 
of the bands near Fermi surface based on the symmetry 
argument in sec. IIIII Furthermore, our model Hamil- 
tonian for the conduction and valence bands is derived 
from the theory of invariants. In sec. IIV( we re-derive 
our model Hamiltonian from the k • p theory and de- 
termine its parameters by more fundamental matrix el- 
ements of the momentum operator in the k • p theory. 
As an application of our model Hamiltonian, the sur- 
face state Hamiltonian and the Landau levels for both 
the bulk and surface states are calculated in Sec. |V] and 
Sec. I VI I respectively. In Sec. IVIH the quantitative limi- 
tation of our model Hamiltonian with four-bands is dis- 
cussed and a new model Hamiltonian is proposed to de- 
scribe the Bi2Ses type of materials more quantitatively. 
In Sec. I Villi we provide a brief discussion and conclusion. 



II. CRYSTAL STRUCTURE, ATOMIC 
ORBITALS AND SYMMETRY 



In this section we will describe the crystal structure 
of Bi2Ses family of materials, and discuss the relevant 
atomic orbit als and the discrete symmetries. A large por- 
tion of the content of this section is already discussed in 
Ref. jl| , but we feel that it is helpful to present the more 
complete version of this discussion here, to make this 
paper self-contained. The crystal structure of Bi2Ses 
is rhombohedral with the space group D^^ {R3m). As 
shown in Fig. [T] (a), the crystal has layered structure 
stacked along z-direction with five atoms (two Bi atoms 
and three Se atoms) in one unit cell, including two equiv- 
alent Se atoms (5'el and 5'el'), two equivalent Bi atoms 
{Bil and BiV) and one Se atom {Se2) which is in- 
equivalent to the 5'el and SeV atoms. Therefore five 
atomic layers can be viewed as one unit, which is usu- 
ally called a quintuple layer. Each atomic layer forms a 
triangle lattice, which has three possible positions, de- 
noted as A, B and C, as shown in Fig(T] (c). Along the 
z-direction, the triangle layers are stacked in the order 

A-B-C-A-B-C . We note that the primitive 

lattice vector U {i = 1,2,3) is not directed along the z 
direction. For example, in one quintuple layer, the Se2 
atoms occupy the A sites; in the next quintuple layer, 
the Se2 atoms do not occupy the A sites but rather the 
C or B sites. Our coordinate is set as the following: the 
origin point is set at the Se2 site; z direction is set per- 
pendicular to the atomic layer, x direction is taken along 
the binary axis with the two fold rotation symmetry and 
y direction is taken along the bisectrix axis, which is the 
crossing line of the refiection plane and the Se2 atomic 



layer plane. The Brillioun Zone (BZ) of this lattice struc- 
ture is shown in Fig. [1] (b). This crystal structure has 
the following discrete symmetries: 

• Three fold rotation symmetry Rs along z direction. 
Rs can be generated by the following transforma- 
tion: X ^ X cos 6 — ysinO, y -> xsinO -\-y cos 6 and 
z ^ z^ where = ^. 

• Two fold rotation symmetry R2 along x direction. 
R2 corresponds to the following transformation: 



Se2 -^ Se2, Bil -^ Bil' Sel 



Sel' 



X ^ x^ y ^ —y. For this symmetry operation, we 
find that Bil (Sel) and BiV {SeV) layers inter- 
change their positions. 

• Inversion P. P: Se2 -^ Se2, Bil -^ Bil' Sel -^ 



Sel' 



-z, X -^ —X, y -^ —y. The Se2 site is 



the inversion center of this lattice structure, hence 
we set Se2 as the origin point. Under inversion 
operation, Bil (Sel) is changed to Bil' {Sel'). 

• Time reversal symmetry T. Time reversal opera- 
tion is given by T = QK, where O = io-2 and K 
is the complex conjugate operator. Here cri,2,3 are 
the Pauli matrice for spin. 




FIG. 1: (a) The crystal structure of Bi2Se^. ti,2,3 is the 
primitive lattice vector, given by ti — (A/3a/3, 0, c/3), t2 = 
( — A/3a/6, a/2, c/3), ts = ( — A/3a/6, —a/2, c/3); where a is the 
lattice constant in the x-y plane and c is the lattice constant 
along the z direction. The quintuple layer is shown in the red 
box with Sel - Bil - Se2 - Bil' - Sel'. Sel (Bil) and Sel' 
{Bil') are equivalent, (b) Brillioun Zone of Bi2Ses. (c) The 
in-plane triangle lattice has three possible position A, B and 
C. 

In order to get a physical picture of the band structure 
of Bi2Ses, we start from the atomic orbitals of Bi and ^e. 
The electron configuration of Bi is 6s'^6p^ and that of Se 
is 4s^4p^. The outmost shells for both Bi and ^e are p- 
orbital, therefore it is natural to consider only p-orbitals 



of Bi and Se and neglect other orbitals. As discussed 
above, Bi2Ses has a layered structure. The chemical 
bonding is very strong within one quintuple layer but 
the two neighboring quintuple layers are only coupled 
by the van der Waals force. Therefore it is reasonable 
for us to first focus on one quintuple layer. Within one 
quintuple layer there are 5 atoms in one unit cell and each 
atom has three orbitals {px^ Py and pz)^ therefore totally 
there are 15 orbitals. The spin is neglected first and 
will be discussed later when we introduce SOC into the 
system. We denote these orbitals as | A, a) with A = Bil^ 
BiV ^ S'el, 5'e2, S'el' and a = p^^ Py^ Pz- As shown in 
Fig(T](a), the Se2 atomic layer stays in the middle of the 
quintuple layer and is sandwiched by two Bi layers {Bil 
and BiV), while two Se layers (5'el and 5'el') are located 
at the outermost. Since all the Se layers are seperated 
by Bi layers, the strongest coupling in this system is the 
coupling between Bi layers and Se layers. Such coupling 
causes level repulsion, so that the Bi energy levels are 
pushed up and form new hybridized states \Ba) and \B^^) 
while the ^e energy levels are pushed down and yield 
three states \Sa)^ \S'^) and |S'Oa), as shown in Fig|2](I). 
Since the system has inversion symmetry, it is convenient 
to combine these orbitals to form the bonding and anti- 
bonding states with the definite parity, which are given 

by 



pi+,±§) 




(I) 



(II) 



(III) 



(IV) 



FIG. 2: Schematic picture of the origin of the band structure 
of Bi2Se3. Starting from the atomic orbitals of Bi and Se, 
the following four steps are required to understand the band 
structure: (I) the hybridization of Bi orbitals and Se orbitals, 
(II) the formation of the bonding and antibonding states due 
to the inversion symmetry, (III) the crystal field splitting and 
(IV) the influence of the SOC. 






mentum: 



(1) 



with the upper index denoting the parity and a = 
Px-,Py,Pz' When the coupling between \Bo^{Sa)) and 
\B'^{S'^)) is taken into account, the bonding and anti- 
bonding states are split, with the anti-bonding state hav- 
ing higher energy than the bonding state. Therefore as 
shown in Figl2](II), the states |Pl+,a) and \P2~ ^a) are 
found to be near the Fermi surface, hence we focus on 
|Pl+,Qf) and |P2~,a) {a = Px^Py^Pz) and neglect the 
other states. Furthermore the crystal has layered struc- 
ture, so the z direction is different from the x or y di- 
rections in the atomic plane. Thus there is an energy 
splitting between Pz and Px,y orbitals for both P1+ and 
P2~ states. We find that \P1~^ ^Px^y) orbitals have higher 
energy than \Pl~^^pz)^ while \P2~ ,px^y) orbitals have 
lower energy than \P2~ ,pz). Consequently, the conduc- 
tion band mainly consists of \Pl~^,Pz) while the valence 
band is dominated by the \P2~ ^pz) orbital before SOC 
is considered, as shown in Fig |2] (III). 

Next we include SOC effect in the above atomic pic- 
ture. The states |Pl+,a,cr) and |P2~,a,cr) are all dou- 
ble degenerate, with one more index a =t,i to de- 
note spin. The atomic SOC Hamiltonian is given by 



H. 



As • L with A 



1 dU 

27tIqc2 r dr 



-^ depending on the de- 



tail potential U of atoms, which couples orbital angular 
momentum to spin. It is convenient to transform the Px 
and Py orbitals to p± with definite orbital angular mo- 



|A,P+,^) 



'V2 



(|A,p^,cr)+z|A,p^,a)), (2) 



|A,P-,cr) = — (|A,p^,cr) -i\A,py,a)), (3) 

where A = P1+, P2~. Within this basis, the atomic SOC 
Hamiltonian is given by 

(A,p+,t \H,o\A,P+,t) = (A,p-4 \H,o\A,P-,i) = ^ 

{A,p+,i \H,,\A,p+,i) = (A,p_,t \H,,\A,p_,t) = -^ 

Aa 



{A,p+,i \H,o\A,p„t) = (A,p_,t \Hso\A,p„i) = — = 



{A,p,,t{l)\H,o\A,p,,t{i)) = 0. 



(4) 



Here the value of A a is a linear combination of the SOC 
coefficient for Bi and Se, depending on how much the 
orbitals of Bi and Se are mixed into the state |A). The 
sign of Aa is always positive for A = P1+, P2~ since the 
potential is always attractive for atoms. As we see, since 
the total angular momentum along the z direction is still 
conserved, hybridization only occurs between lA,^;^,^) 
{\A,pz,i)) and |A,p+,|) (|A,p_,t)). After taking into 
account SOC, the new eigen-states are given by 



|A,-) = |A,p+,t) 



|A,- 



|A,p-,;) 



(5) 
(6) 



\A+,-) =u^\A,p,,t) + vi\A,p+,i) 
|A_,i) = u^|A,p„t> + «^|A,p+,;) 



(7) 
(8) 



|A+,--) = «)*|A,p„;) + («^)*|A,p_,t> (9) 

|A_,-i) = («^r|A,p„;> + («^r|A,P-,t> (10) 

with the eigen-energies E^,^ and E^^ (each is double 
degenerate) and u, v obtained by solving the following 
2x2 Hamiltonian 



ff^f Ek,x - Aa/2 Aa/\/2 \ 
\k/V2 Ea,, J 



(11) 



For the above eigen-states (j5j)'^([T0]), all the information 
about SOC is included in the coefficient u and v, which 
are given by 



1 



A^A± 



(A^a)' + ^ 



Aa/V2 



(12) 



A2 



2AEl 



± 



Ea,x — Ea,z—>^a/'^ 



explicitly, where N± = 

2AEa^AEI^XI/2 and A^;^ 
The energy splitting between the Px{y) orbital and 
the pz orbital due to the crystal field is larger than 
the energy scale of SOC and A£^a is dominated by 
Ea,x — Ea^z' Now as we see, the SOC couples lA,^^,^) 
{\A,Pz,i)) to |A,p+,|) (|A,p_,t)) so that it induces the 
level repulsion between these two states. Consequently, 
|Pll,±|) is pushed down while |P2^,±^) is pushed 
up, which yields the level crossing between these two 
pairs of states, when the SOC is strong enough, as shown 
in Fig 12] (IV). Since these two pairs of states have the 
opposite parity, their crossing leads to a band inversion^ 
similar to the case in the HgTe quantum wells [5]. This 
is the key signature of the topological insulator phase in 
Bi2Ses family of materials [1]. Therefore in the following 
we will focus on these four states and regard the other 
states as the perturbation. 



III. MODEL HAMILTONIAN DERIVED FROM 
SYMMETRY PRINCIPLES 

From the discussion of the atomic orbitals in the last 
section, we obtain an intuitive physical picture of the 
band structure of Bi2Ses. Compared with the ab initio 
calculation, we can denote the bands near Fermi surface 
by |A^,a) where A = P1±^P2± and a = ±^,±|, as 
shown in Fig [3l Roughly, these states mainly consist of 
the bonding or anti-bonding states of the p-orbitals of Bi 
or Se atoms. However, other orbitals such as s-orbitals 
of Bi and Se will also mix into these states. To iden- 
tify each band without any ambiguity, it is necessary to 
relate each band with the representation of the crystal 
symmetry. At F point, each state should belong to an 



irreducible representation of the crystal symmetry group 
and the hybridization between orbitals preserve the sym- 
metry properties. Therefore, a suitable method to iden- 
tify each band is to use the symmetry of the crystal. In 
this section, we will first identify each band according to 
the irreducible representation of the crystal group D^^ 
and then try to derive our model Hamiltonian just from 
symmetry principles. 

First let's consider the states without spin, which are 
denoted as |A^,a) with A = P1,P2 and a = Px^Py^Pz- 
The crystal of Bi2Se^ belongs to the group D^^ with 
the character table given in table @ of Appendix [A[39| . 
Since the crystal is inversion symmetric, each represen- 
tation has a definite parity eigenvalue. For each parity, 
there are two one-dimensional representations F^ and 
V2 and one two-dimensional representation Fg , where 
the upper index denotes the parity (+ for even and - 
for odd). According to the wave functions constructed 
from the simple atomic orbital picture, we can deter- 
mine the transformation property of the wave functions 
under the generators i?3, R2 and P of the point group. 
For example, let's look at the operation R2 on the state 
|Pl+,j9a^) = -^{\Bx) - \B'^))- The R2 rotation does 
not change the px orbital, however it changes the po- 
sition of Bil {Sel) and BiV {SeV) and correspondingly 
changes \B) to \B')^ thus we should have R2\P1^ ^Px) = 
— |Pl+,Pa;). Similiar argument can be applied to other 
states and finally the transformation of the states under 
the cyrstal symmetry operation is listed as follows. 

• Three fold rotation symmetry Rs: \A^^px) -^ 
cose\A^,Px) - s'mO\A^,py), \A^,Py) ^ 
smO\A^,Px) + cosO\A^, py) and lA^,^^) -^ 
|A±,p,),with^=f. 

• Two fold rotation symmetry R2: \A^,px) -^ 
T|A±,p^), \A^,Py) -^ ±|A±,py), \A^,Pz) -^ 
±|A±,p.). 

• Inversion P: \A^ ^a) ^ ^lA"^^ a)^ a = Px^Py^Pz- 

Here A = P1±,P2±. According to the above transfor- 
mation, we find that \A'^'^~\px) and |A+*^~\p^) belong 
to the fg ^"^representation. \A~^^pz) belongs to f^ rep- 
resentation and lA",^^) belongs to F^ representation. 

To take into account spin, we introduce the spinor rep- 
resentation Fg", which changes its sign under the rota- 
tion C = 27r. The double group of D^^ can be con- 
structed by the direct product of f f 2 3 ^^^ ^g". As 
shown in (|AT]) ^ (|A3|) . we find that f ^ (g) f^ will give 
two new one-dimensional representations 1^4 and T^ , 
which are conjugate to each other. The character ta- 
ble of the double group for D^^ is given in table (jTT]) 
of Appendix [A[39]. With SOC, the eigen-states in (|5j) 
~ (|6]) can also be analysed by the decomposition of di- 
rect production. From (|A2[) and (|A3[) . the direct product 
of Fg and r^2 always gives Fg representation, therefore 
|A+, ±^) with A = PI, P2 should belong to f^ represen- 
taion while |A~, ±^) should belong to F^ representation. 



The states |A^,zb3/2) originate from the combination of 
\K^Px,y) and spin. According to (|Aip , it is expected that 
|A^,±3/2) should be a combination of fj and f^ rep- 
resentations. Indeed by carefuhy inspecting the transfor- 
mation behavior under the operation R2 and i?3, we find 
that 



|A±,f4) = -i=(|A±,3/2) + |A±,-3/2)) 

belongs to the F^ representation, while 

|A±,f5) = -^(|A±,3/2)-|A±,-3/2)) 



(13) 



(14) 




l^i|.±5) r^ 



P2-,±|)r4r5 
P2i,±i)r7 



belongs to the f^ representation. The above results can 
also be worked out by considering the forms of the trans- 
formation for the states (|5]) ~ (|6]), which are given by 



Three fold rotation symmetry R^: |A, ±^) 



-^ 



e±^t|A,±i), |A,±f) 
P1±,P2±. 



|A,±|), where A 



Two fold rotation symmetry R2: |A+,±^) -^ 
z|A+,T^), |A-,±^) ^ -^|A-,T^), |A+,±f) ^ 

i|A+,Tf), |A-,±i) - 
P1±,P2±. 



-z|A",^f), with A = 



• Inversion P: |A^,a) -^ ±|A^,a), with A = 
P1±,P2± anda = ±|,±^. 

It is instructive to compare the present case with the 
more common semiconductor crystal structures, such as 
diamond or zinc-blende structure. In that case, the cou- 
pling between p orbitals and the spin usually gives the 
four-dimensional Fg and the two-dimensional F7 repre- 
sentation. In the present case, due to the lower symme- 
try of the crystal structure, the Ft representation is the 
same as the Tq representation while the Fg representation 
is reduced to two one-dimensional representations r4 and 
F5 and one two-dimensional representation Tq. In Fig[3l 
the representation of the bands near the Fermi surface is 
given. 

Next we derive our model Hamiltonian to describe the 
low energy physics of Bi2Ses just based on the symmetry 
of the wave function at F point. As described above, 
near the Fermi surface the conduction band and valence 
band are determined by the four states |Pll,±^) and 
|P2^,±|), belonging to the Tq and F^ representation. 
Therefore the minimum model Hamiltonian for Bi2Ses 
should be written with these four states as the basis. 
Generally, any 4x4 Hamiltonian can be expanded with 
Dirac F matrices as 



Heff = e(k)I + ^ di{k)ri + J2 dijO<-)n 



(15) 



where I is the 4x4 identity matrix, F^ (z = 1 • • • 5) de- 
note the five Dirac F matrices satisfying {F^,Fj} = 2(5^j, 



FIG. 3: The band structure of Bi2Se3 is obtained from ab 
initio calculation, and the bands near Fermi surface is iden- 



tified with |A^,a). A = P1±,P2± and a — ^^,^^. 
corresponding irreducible representation is also given. 



±i 



The 



and the ten anti-commutators of F matrices are given 
by Tij = [Fi,Fj]/2z. e(k), (ii(k) and dij{k.) can be ex- 
panded by the powers of the momentum k. The con- 
struction of F is given in Appendix [Bl Now let's assume 
the above Hamiltonian is written in the basis |Pll, ^), 
|P2^,|), |Plt,-^) and |P2^,-|). Then according to 
the transformation of the states under the symmetry op- 
eration discussed above, we can construct the transfor- 
mation matrices as follows. 

• Time reversal symmetry: T = QK^ where Q = 
ia2 1 and K is the complex conjugate operator. 

• Three fold rotation symmetry operation: Rs = e*^^ 
with H = 0-3 1 and 6> = ^. 

• Two fold rotation symmetry operation: R2 = i(Ji 
T3. 

• Inversion: P = 1 r3 . 

In the above, cr acts in the spin basis and r acts in the ba- 
sis of P1+ and P2~ sub-bands. According to the above 
transformation matrices, we can obtain the irreducible 
representation of each F matrix, details of which are de- 
rived in Appendix [Bl The invariance of the Hamiltonian 
requires that the function (ii(k) {dij{k)) should have the 
same behavior to the corresponding F^ {Tij ) matrix under 
the symmetry operation, which means that they should 
belong to the same representation of the crystal point 
group. In the table [Till of the Appendix [Bl we list the 
representation for both the F matrices and the polyno- 
mials of k, and also their transformation properties under 
the time reversal operation. Since we hope to preserve 
both time reversal symmetry and crystal symmetry, we 
must choose the F matrices and polynomials of k with 
the same representation. For example, Fi and F2 carry 
the representation F^ and are odd under time reversal. 



and so are kx and ky. Therefore they can together form 
an invariant term for Hamiltonian. Finahy, up to 0(/c^), 
our model Hamiltonian yields 



^eff = Ho 



H', 



H'o 



e^,^M{\^)V^^B{k,)V^k, 
-^A(k\\){Tiky-T2kx) 



(16) 



H'^ = RiTsikl - Skxkl) + R2T43klky - k^) (17) 

where ek = Co + Cife^ + C2A;|, and A1(k) = Mo + Mi/c^ + 
M2kj and^(A;||) = Ao^A2kj and B{kz) = 5o+^2^^ and 
k'u = kl ^ ky. Hq keeps the in-plane rotation symmetry 
along z direction, while H^ breaks the in-plane rotation 
symmetry down to three fold rotation symmetry. In the 
following, we find the bulk H^ term will also lead to the 
correction to the effective surface Hamiltonian, which has 
been studied in Ref. [35]. 

The above Hamiltonian is the same to that presented 
by Zhang et al in the Ref. [l|, which can be shown by 
performing the transformation 



Ui 



10 
-i 
10 
i 



(18) 



and the Hamiltonian is transformed into 
Heff = Ho + H3 
Ho = UiH'oUl = ek + 

/ M{\i) B{k,)k, 

B{k,)k, -7W(k) A{k\\)k- 

A{k\\)k+ M{\i) 

\^(fc||)fc+ -B{k,)k, 



H3 = UiH^ul 



R2{kl - fc3 ) 



Ri{kl + k^_) 



-i 
-i 
i 
i 




(20) 



Now we can see that Hq is nearly the same to the model 
Hamiltonian (1) in Ref. [l|, except the A2 term and B2 
term, which represent the high order correction to the 
Fermi velocity Aq and Bq. Since such correction is not 
important near F point, we will neglect these two terms 
in the following. Derivation of our model Hamiltonian 
([19]) and ([2Q|) from the symmetry principles is the central 
result of this section. 



IV. MODEL HAMILTONIAN DERIVED FROM 
THE k p PERTURBATION THEORY 

Up to now we have obtained our model Hamilto- 
nian from the symmetry principles, or the theory of 



invariants [40]. In this section, we will derive the model 
Hamiltonian through another way, k • p theory, and con- 
nect the parameters of the model Hamiltonian to the 
more fundamental matrix elements of momentum in k • p 
theory. 

The basic idea of k-p theory is to use the wave function 
at F point in BZ as the zeroth-order wave function and 
treat 



H' 



h 

mo 



kp 



(21) 



as a perturbation, where p = —ihdr is the momen- 
tum operator acting on the zeroth-order wave function 
and the crystal momentum k is regarded as a small pa- 
rameter during the perturbation procedure. The model 
Hamiltonian is expanded by the powers of k. With the 
perturbation formalism (jClSp ~ (jC2ip . we can project 
the system into the subspace spanned by the four states 
|Pli,l/2) ^ |1), |P2;,l/2) ^ |2), |Plt,-l/2) ^ |3) 
and |P2^,— 1/2) = |4), which are used as the basis of 
our model Hamiltonian. All the other states are treated 
in the perturbation procedure and the details are given 
in Appendix [O The obtained model Hamiltonian will 
depend on a series of matrix elements of momentum 
(Al, a|p|A2, /3), which can be simplified due the symme- 
try of the crystal. For example, due to the existence of 
the inversion symmetry, all the states at F point have 
definite parity eigenvalues. Since the momentum p has 
odd parity, the matrix elements of momentum between 
two states with the same parity always vanish. The wave 
function at the F point can be obtained through ab inito 
calculation, consequently all these matrix elements can 
be calculated. With these matrix elements, we apply the 
perturbation formulism (jClSp ~ ()C2ip to the system and 
recover our model Hamiltonian ([19]) and (|2Q|) . The pa- 
rameters of our model Hamiltonian Co, Ci, C2, Mq, Mi, 
M2, Ao, Bo, Ri and R2 can be expressed as the function 
of the parameters Pa A,, Q a^^ a^, M k^m^ ^Ai,A2^ ^Ai,A2 
and 5'ai,A2 through ([cS2]) - ([UgT]) in the Appendix O 
With these expressions, we can numerically calculate the 
values of the parameters of our model Hamiltonian, which 
is listed in table llVi We note that the parameters given 
here are different from those in Ref. [1] and [41], where 
the parameters are determined by fitting to the energy 
dispersion, which has some ambiguity. In the present 
method, since we directly calculate the matrix elements 
of momentum from microscopic wave functions, there is 
no ambiguity. 

The key result of this section are the parameters of our 
model Hamiltonian, given in table IIVI The fitted energy 
dispersions for Bi2Ses, Bi2Tes and Sb2Tes are plotted 
in FigHl It is shown that our model Hamiltonian is valid 
in the regime kx,z < 0.04A~^. However the maximum of 
the valence band for Bi2Ses and Bi2Tes stays away from 
F point, at about kx ~ 0.07A~^, therefore we need to 
keep in mind that there may be some discrepancies when 
we try to use our model Hamiltonian to describe Bi2Ses 
quantitatively. Since our model Hamiltonian with four 
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FIG. 4: The energy dispersion obtained from our model 
Hamiltonian with four bands (sohd hne) is compared with 
that from ab initio calculation (dashed line). Here (a) and 
(b) is for Bi2Se3, (c) and (d) is for Bi2Tes, while (e) and 
(f) is for Sb2Te3. In (a), (c) and (e), the red line represents 
the dispersion along kx direction while the blue line is for ky 
direction. 



bands already captures the salient features of the band 
dispersion, especially the inverted band structure, in the 
following two sections, we still stay in the framework of 
our model Hamiltonian to discuss the topological surface 
states and the Landau levels in the magnetic field. Then 
in the last section, we will extend our model Hamilto- 
nian to include eight bands, in order to describe these 
materials more quantitatively. 



V. SURFACE STATES 

An important physical consequence of the non-trivial 
topology is the existence of topological surface states. In 
this section, we would like to study the surface state and 
its effective Hamiltonian based on our model Hamiltonian 
derived above with open boundary condition. 

Consider our model Hamiltonian ([T7j) defined on the 
half space given hy z > 0. We can divide our model 



Hamiltonian into two parts 

H = Ho^Hi (22) 

Ho = i{h) + M{k,)r^ + BoT^k, (23) 

-^Ao{Tiky-T2k,)^Hs. (24) 

where e{kz) = Co + Cik^ and M{kz) = Mq^ Mik^. All 
kz dependent terms are included in Hq. We replace kz 
by —idz, and obtain the eigen- value equation 



Ho{kz^-idz)^{z) = E^{z) 



(25) 



Since r4 = 10T2 and Fs = l^rs are both block-diagonal, 
the Hamiltonian Hq is also block-diagonal and the eigen- 
states have the form 



^t(^) 







^iiz) 





^0 J' 



(26) 



where is a two-component zero vector, ^^{z) is related 
to "^liz) by time-reversal operation. To obtain the sur- 
face states, the wave function ^|Jo{z) should be localized 
at the surface and satisfies the eigen-equation 

{i{-idz) + M{-idz)rs - iBoT2dz) M^) = ^V^o(^)(27) 

which has been solved analytically for the open boundary 
condition with different methods. [31, 32, 41-43] Here in 
order to show the existence of the surface states and to 
find the region where the surface states exist, we would 
like to breifly review the derivation for the explicit form 
of the surface states by neglecting the e for simplicity [3l|. 
After neglecting the e term, the eigen equation ([27|) 
has the particle hole symmetry, therefore we would like 
to expect a sepecial surface states with E = can exist. 
With the wave function ansatz ipo = (j)e^^ ^ the above 
equation can be simplified as 

(Mo - MiA^) n^ = Bq\(I). (28) 

It is obvious that the two-component wave function (j) 
should be the eigen-state of Pauli matrix ri. Let's define 
Ti(j)± = ±^±, then the equation ([28|) is simplified to a 
quadratic equation for A. Another important observation 
is that if A is a solution for ^+, then —A is the solution 
for (j)- . Consequently, the generic wave function is given 

by 



ipo{z) = {ae 



\\Z 



he 



Mz\ 



+ (ce" 



-\\Z 



-de-^2^)(^_,(29) 



where Ai^2 satisfy 



Al,2 



^(-BoivWwTfB^y 



(30) 



Similar to the Ref. [31], the open boundary condition 
?/^(0) = 0, together with the normalizability of the wave 



function in the region z > 0, leads to the existence con- 
dition of the surface states, 3?Ai^2 < (c = d = 0) or 
J?Ai,2 > (a = 6 = 0), which can only be satisfied with 
the band inversion condition MqMi < 0. Furthermore, it 
is easy to show that when Bq/Mi > 0, 3?Ai,2 < 0, while 
Bo /Ml < 0, 3? A 1,2 > 0, thus the wave function for the 
surface states at T point is given by 



i^o{m) 



a (e 
c (e 



\\Z 



(^^-\iz 









5o/Mi < 



We emphasize here the sign change of Bq/Mi will change 
the spin basis of the surface states, which will be the 
key point to determine the helicity of the Dirac Hamilto- 
nian for the topological surface states. Another impor- 
tant quantity of the surface states is the decaying length, 

which can be defined as Ic = max \ ^otv — \t \ , given by 



lc= { 



2 Mo 



Bo-y^4MoMi+Bg 
2 Mo 



Bo + v'4MoMi+Bg 
2Mo 




(32) 



Bo-y^4MoMi+Bg 
2Mo 



where 5R takes the real part. 

In the above, we take a simple tight-binding model to 
show the existence condition and the form of the wave 
function for the surface states, which can help us to 
understand the underlining physics qualitively, but not 
quantitively. In the realistic materials, the detail form of 
V^o will depend on the material detail, such as the bound- 
ary condition or the detail parameters, however the form 
of the wave function (|26]) remains valid. Therefore in the 
following, we just simply treat i/jq by some parameters. 
On the sub-space ^ = [^^, ^|], we find that 



(*|ri|*> 


= aic^x, 


{*|r2|*> = 


= a^a. 


(^irsl*) 


= aiCTz, 


{*|r4|*) = 


= 0, 


(^irsl*) 


= "3, 







(33) 



with ai = (V^olnl^o) and as = (V^okalV^o)- With 
these expressions, the effective Hamiltonian of the sur- 
face states ^ is given by[l| 



Cok^ 



Hsur — Co -f- L/2 

= (Co + agMo) 



A{aa:ky - (jyk:,) + R(k% + k^_)(j^ 

(C2 + asM2)k\ + A^aiiycrxky - CFykx) 

kl)cj, (34) 



with k± = kx ^ iky = /cye^*^. The k^ terms have also 
been found in Ref. [35| . In the following numerical calcu- 
lation, the coefficient ai and 0^3 are treated as two fitting 



parameters to the experiment, given by ai 

Cexp~Co 



0.99 



and as 



Mo 



-0.15, where Ae 



A 

-- 3.29eV • A 



comes from the Fermi velocity of the surfaces states and 
Cexp = 0.035ey comes from the position of the surface 
Dirac points [10|- Further we need to check the spin op- 
erators in this system. Again we use the wave function 
from ab initio calculation and project the spin operator 
into the subspace spanned by the four basis states. After 
we obtain the spin operators for our model Hamiltonian, 
we can use the eigen wave function (|26|) to project the 
spin operator into the surface states subspace. Finally 
we find that (^l^^^lvE^) = S^octx, {"^ISyl"^) = SyQCjy and 
{^\Sz\^) = SzoCTz with Sx{y^z)o to be some constants. 
This indicates that a matrix in our model Hamiltonian 
([34|) is proportional to the real spin. 

The derivation of the surface Hamiltonian (|34|) is the 
central result of this section. In the limit /c ^ 0, the 
linear term in Hamiltonian (|34|) will be dominant, then 
the surface states show the linear dispersion with heli- 
cal spin texture, which has the opposite direction for the 
conduction and valence band, as shown in Fig|5](a). Such 
type of spin texture is similar to one of the fermi surfaces 
in the usual 2D electron gas with Rashba SOC[40|, 1441 1. 
which can be simply understood from the fact that the 
inversion symmetry is broken near the surface. The he- 
lical spin texture has also been calculated from ab initio 
method [45] and already observed in the pioneering spin- 
resolved APRES measurement [18]. From (|34|) . the helic- 
ity of the spin texture is determined by the sign of the 
coefficient A = A^ai^ where ai = (V^okil'^o) is related to 
the spin basis of the surface states (|3T]) . Therefore the he- 
licity is determined by the relative sign of Aq and Bq/Mi. 
Furthermore, due to the inversion condition MiMq < 0, 
the sign of Mi is already determined by the gap Mq. 
Consequently, within our model Hamiltonian the helicity 
of the spin texture is given by the relative sign of the 
coefficients of two types of linear terms, Aq and Bq. 

To further explore the origin of the helicity of the spin 
texture in the atomic levels, we relate the coefficients 
Ao and Bq to the atomic SOC by using the expression 
dlD-^dlOl), as 



2mn 



(«^i\f")*(Pl+,p,b+|P2-,p_) 



+ (v^i^uf ")* (Pl+,p+b+|P2-,p,)' 

Bo = —{Plt,Pz\Pz\P2-,p,) 
mo 

= — \{u^'^ru^^~{Pl+,Pz\Pz\p2-,Pz) 
mo L 

+(t;^i^)*t;f"(Pl+,p+b.|P2-,p+)". 



(35) 



(36) 



Here |A,a) (A = P1+,P2~ and a = Px^Py^Pz) are the 
atomic orbit als without any SOC and all the dependence 
of SOC are included in the coefficient u^ and v^. From 
(p!2|) . we find that for u^^v^^ it is proportional to Aa^ (or 
AA2), which indicates that Aq depends on the sign of the 



atomic SOC, while it is only possible for {u^^Yu^'^ and 



{v 



Ai\*.,A2 



to be independent of Aa^ 2 ^^ be proportional 



to Aa^Aa2, thus the sign of ^o will not depend on the 
atomic SOC. Finally, we conclude that the helicity of the 
spin texture is originally related to the atomic SOC. 

In the above, we have shown how the the linear term 
determines the spin texture of the surface states, which 
can also be affected by the quadratic term and cubic term 
in the effective Hamiltonian (|34|) . We solve the eigenvalue 
problems of the whole effective Hamiltonian and obtain 
the eigen-energy and eigen states as 



E± = Co^ C2kl ± \ A^kl + m^k^ cos2 3(9 (37) 



^± = 



Vat \d±-R{k\ + k^_) 



(38) 



with d-\ 



±Ji2A:2+4^2^6cos2 3(9 and TV = A^k'^ 



( \/i2p^4^2^6cos2 3l9 - 2Rk^ cos 3(9) 2 . Consequently 
the spin polarization in k space is given by 






2Ak 
2At 



^((i+ -2Rk^cos30) (39) 

(d+ -2Rk^cos30) (40) 



N 



(V^+k.|V^+) = 1^^!^(^^ - 2m^cos3^)(41) 

which is plotted in Fig[5](b). In the limit /c ^ 0, the spin 
polarization almost lies in the xy plane, which is due 
to the linear term and has been discussed in the above. 
When k is increased, the k-cubic term comes into play, 
which will not only induce the hexagonal warping of the 
constant energy contours [35| but also yield z direction 
spin polarization, similar to the situation in Bi^Sbi-x 
studied by ab initio calculation [46.]. 




(b) 




Topological insulators 



FIG. 5: (a) Spin texture of the surface states near F point. 
For conduction band, the helicity is left handed while for va- 
lence band, it is right handed, (b) Spin texture of the con- 
duction band of the surface states in the momentum space. 
The arrow represents the x-y planar spin polarization while 
the color indicates the z component of the spin polarization. 
Here red is for spin up while blue is for spin down. The black 
line gives the constant energy contours. 



> 

LU 




2 4 6 8 

sgn{N)^NB (T5) 



VI. MAGNETIC FIELD AND LANDAU LEVEL 

In this section, we study the Landau level problem for 
both the bulk states and the surface states, which is im- 
portant for predicting or understanding many properties 
of the system in a magnetic field, such as the SdH oscil- 
lation, surface quantum Hall effect and magneto-optics. 
In this regards, our model Hamiltonian has the unique 
advantage, since the magnetic field effect cannot be in- 
corporated in ah initio calculations. For the realistic fi- 
nite sample, the bulk Landau levels will always coexist 
with the surface Landau levels, thus both the two types of 
Landau levels need to be taken into account. For simplic- 
ity, we solve the bulk Landau level for an infinite sample 
and the surface Landau level for the semi-infinite sam- 
ple. The mixing effect between bulk and surface Landau 
levels is neglected here. 

For the bulk state, there are two types of contribution 
from the magnetic field, the orbital effect and the Zee- 
man effect. The orbital effect can be included by Peierls 



FIG. 6: The energy of the Landau level verse sgn{N)\NB 
where N is the Landau level index and B is the magnetic 
field. 



substitution [475 k ^ tt = k+|A with A = (0, B^x, 0) for 
magnetic field along z-direction. We introduce the anni- 
hilation and creator operators a 



with L 



for the harmonic oscillator function (f^i 



a and a' satify a(^N 



N^. 



N-l 



a^ifN 



VN + IfN+l 



and [a,a^ = 1. With the operators a and a^ , the Hamil- 



tonian (1191) is written as 



H< 



OB 



( M Bok, 

Bok, -M 

Ao^at 
\Ao^a^ 



Ao^a\ 



M 
-Bok, 



la 


-Bokz 
-M 



(42) 
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FIG. 7: The Landau levels in the magnetic field for both 
the bulk states (red lines) and the surface states (blue lines) 
are shown in (b). Also the dispersion along z direction at (a) 
B = 0T and (c) B = 20T is plotted. 



= Co 



2M2 



where e{kz,a^a) 

M{kz,aU) = Mo^Mikl + ^ (ata + \). The k^ term, 
which breaks the in-plane rotation symmetry, is neglected 
here, therefore the wave function should have the form of 
^AT = [fi ^N-iJ2 ^N-1, f^ ^N Ja ^nV • With such 
wave function ansatz, the Hamiltonian is transformed to 



HoB{h,N) 



with M^ 



Bok, 




Bok 
M 



N-\ 

■2N 







Aa 



iO- 



2N 







-Bokz 



/2N 


-Bok^ 



;43) 



/ 



e{k,,N) 



e{kz,N) and Mp^ 
i{k.,N). 

To consider the Zeeman splitting, we need to further 
calculate the effective g-factor|4Qll48|. The effective Zee- 
man type coupling can also written down for our model 
Hamiltonian by symmetry principles. By a quick inspec- 
tion of the table piip , we find that the following terms 



Hz = gzl^l2Bz + gz2^MBz 



-^gxy2 ( ri4 T2A 



COS (pi sm 01 

— sin (/>! cos (/>! 

cos (j)2 sin 02 

— sin 02 cos 02 



Bo, 



bI J (44) 



are possible couplings to the magnetic field. Again 0i and 
02 are phase factors which need to be determined from 
other methods, and here we take 0i = 02 = to coincide 
with the results from the k-p method, and explicitly (|i^ 
can be written as 



gizBz gipB_ 

Mb I g2zBz g2pB- 

2 I gipB^ -gizBz 

g2pB^ -g2zBz 



(45) 



with i^B = 27^, and g^yi + gxy2 = ^gip, gxyi - gxy2 = 

^g2p, gzi^gz2 = ^giz, gzi-gz2 = ^g2z- This model 



Hamiltonian can also be derived from k • p theory and 
the parameters ^i^, g2z-, gip and g2p can be related to 
the matrix elements of the momentum operator p in the 
k • p theory, with full details given in Appendix O Now 
our total Hamiltonian for the bulk states under the z- 
direction magnetic field is given by Hb = Hqb + Hz, 
which can be solved numerically to obtain the Landau 
level E^^^^{B^kz) with Landau level index N and band 
index rj under the z-direction magnetic field. 

Similar procedure can be applied to the surface effec- 
tive Hamiltonian (|34|) . With the wave function ansatz 
^ sur.N = [gi ^N-i,g2 ^nV , the surfacc Hamiltonian is 
changed to 



J^sur,BO[^^ ) = <-^0 H T O2 r (T^ 



2NeB^~ 



(46) 



and and the Zeeman type term is given by 



llb Ub 

Hsur,Z = —gszO'zBz + -^gsp {^xBx + ^yBy) (47) 

with ^gsz = gzi^gz2(^2. and ^gsp = gxyi^gxy2(^3' The 
total Hamiltonian for the surface states yields Hsur,B = 
Hsur,BO + Hsur,z and correspondingly the Landau level 
in z-direction magnetic field Bz is solved as 



^r^(7V) = Co + ^^C2- 



'fi- 



C2eBz I /£B 

h ^ 2 



gszBzj 



2NeB, 



A^ (48) 



with 5 = ± for A/" = 1, 2, • • • and 



£;.nr(Q) ^ ^^ 



eB^ 



C2 



rBy 



(49) 



for zero mode A/" = 0. Here we note that due to the ex- 
istence of the quadratic term C2 k'i , the square root de- 
pendence of the energy level verse magnetic field is only 
an approximation applicable for low Landau levels and 
low magnetic field. For high magnetic field, it will be a 
combination of the linear contribution and square root 
contribution. As shown in Fig [6l the energy of the Lan- 
dau levels are plotted as the function of sgn{N)\/NBz 
and the non-linear behavior will appear for high ^/WB^. 
In Fig [3 (b), the Landau levels for both bulk states 
and surface states are plotted as a function of magnetic 
field. Here we emphasize that for the bulk states, the 
Landau levels are plotted for the /c^ = point (red line 
in Fig [71 (b)). The dispersion along z direction is also 
shown in Fig [71 (a) for B = OT and (c) for B = 20T. 
From Fig [3 (a), we find that the maximum of the valence 
band is not located at /c^ = point for small magnetic 
field, hence we plot the maximum of valence band as 
the green lines in Fig [71(b). When the magnetic field is 
increased, the bulk gap is also decreased by a significant 
amount (about 160meV for 5T magnetic field), which is 
due to the double hump structure for the valence band 



11 



S 

s 
8 



M 

■i' 




bulk 

^ surface 

total 


\\a/W^ 



-0.2 



0.2 

E(eV) 



0.4 



FIG. 8: The density of states (DOS) as a function of energy is 
plotted with B = lOT for the bulk (blue line), surface (green 
Une) and total (red line) Landau levels. 



> 
3, 

LU 



> 

S 

LU 



(b) 



0.02 0.04 0.06 0.08 0.1 0.01 0.02 0.03 0.04 0.05 

k„(1/A) k(1/A) 



(C) 



> 

S -0.1 



> 

S-0-1 

LU 



(d) 



0.02 0.04 0.06 0.08 0.1 0.01 0.02 0.03 0.04 0.05 



k„ (1/A) 



k ri/A) 



dispersion of Bi2Se^. Such decrease may be observed in 
a magneto-optical measurement. 

In order to compare with the scanning tunneling mi- 
croscope (STM) experiment |22l. |23] . it is helpful to inves- 
tigate the local density of states (LDOS) at the surface. 
The LDOS for the surface states and the bulk states can 
be obtained by|4Q| 
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Dsur{E,B) = 2_^ 
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N,s 



2r2 



27rr 



(50) 



and 



Dbuik{E,B) = Lq —— 



dkz 

2^ 



G 



{E-El 



HB,k,)y 



N,s 



27tT 



-(51) 



respectively, where G = |^ = ^^ is the degeneracy 
of each Landau level and T is the broading. In order 
to compare the bulk LDOS with the surface LDOS, we 
require to introduce a length scale Lq which represents 
the detection depth of STM. Here we simply take Lq 
to be the thickness of one quintuple layer. Furthermore 
the surface states only exist near k = 0, thus we need 
to take a cut-off for the Landau level index N. With 
the formula (|50j) and dSU, LDOS for both the bulk and 
surface Landau levels are shown in Fig. [51 The bulk 
LDOS shows a gap of about 0.3eV and within the bulk 
gap, only surface LDOS remains and shows clearly the 
Landau levels as discreted peaks. The largest Landau gap 
for surface states is between the 0th and 1th Landau level, 
about 50meV, which is large enough for the observation 
of the topological magneto-electric effect |36l, l37J . 



FIG. 9: The energy dispersion obtained from the new model 
Hamiltonian with eight bands (solid line) is compared with 
that from ab initio calculation (dashed line) for (a), (c) and 
(e) kx and ky direction and (b), (d) and (f) kz direction. 
Here (a) and (b) is for Bi2Se3, (c) and (d) is for Bi2Te3, 
while (e) and (f) is for Sb2Te3. In (a), (c) and (e), the red 
line represents the dispersion along kx direction while the blue 



VII. NEW MODEL HAMILTONIAN WITH 
EIGHT BANDS 

As we have described above, our model Hamiltonian 
can capture the salient topological features of the Bi2Ses 
family of materials. However, for full quantitative fitting 
with the first principle calculations, we need to expand 
the basis set. By inspecting carefully the k • p matrix 
elements, we find that there are strong couplings be- 
tween the state |Pll,±|) and the state |P2~,r4^5) or 
|P2I,±|). For example, at the valence band maximum 
kx ^ 0.07A~^, we find that these couplings can be as 
large as the energy gap between these states. Therefore 
it is not surprise that our model Hamiltonian with four 
bands is not suitable in this regime. The strong cou- 
plings between these states indicate that if we want to 
describe this material more accurately, we need to fur- 
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ther include the states |P2 ,r4^5) and |P2_,±^) into 
our model Hamiltonian. In the basis sequence |Pll, ^), 

|Pil,-i), |P2;,i), |Pi;,-^), |P2-,f4), |P2-,f5), 
I 



|P2I, ^) and |P2I, — ^), following the similar perturba- 
tion procedure, we find that our model Hamiltonian is 
written as 
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-iP^k- 

335 (k) 

53*6(-k) 
/5(k) 



hQ2k- 

536 (k) 

-53*5(-k) 



/5(k) 



|fczQ3 

k3*fc+ 

/37(k) 

547 (k) 

-56*8(-k) 

k) 

/7(k) 



jfc-f3 \ 

-54*7(-k) 

/3*7(-k) 
558 (k) 
568 (k) 


/7(k) 



(52) 



with 



VIII. CONCLUSIONS 



/i(ii)(k) 



*(*i) 2 



• ^i(ii)^|| 



^i^(k) = Uijk^k^ + Fi^/c? 



(53) 
(54) 



The parameters Fij , i^^j , Uij and Fij can now also be de- 
termined by the perturbation theory, which is shown in 
appendix O In this Hamiltonian, time-reversal symme- 
try is already satisfied. Furthermore, R2 rotation symme- 
try yields that U35 = ^3*6. ^35 = -^3*6. ^58 = -% and 
V58 = Vq^. The obtained parameters are listed in table 
IVl and the band dispersion is found to fit well with that 
of ab initio calculation, as shown in Fig [9l This demon- 
strate that the eight band model is suitable to serve as 
the basis of the quantitative study of Bi2Ses family of 
materials. 

We would like to make some more remarks about the 
eight band model. Firstly, in our model Hamiltonian 
with four bands, the leading term that breaks the in- 
plane full rotation symmetry down to Rs symmetry is 
the third order in perturbation, while in the eight band 
model, it is the second order coupling gij. These type 
of terms exist because the states |P2~,r4^5) themselves 
break the rotation symmetry according to the expres- 
sion (p!3|) and ([T4j). Secondly, it is interesting to compare 
the present eight band model with the well-known Kane 
model for usual III-V or II- VI group semiconductors with 
zinc-blend structure. In fact there is a one-to-one cor- 
respondence between the basis of these two models, in 
which |Pll,±^) corresponds to the electron band (Fe), 
|P2^,±^) and |P2~,F4^5) correspond to the light hole 
and heavy hole bands (Fg), respectively, and |P2I,±|) 
corresponds to the spin-orbit split-off band (F7). There- 
fore, from the symmetry point of view, our model here is 
nothing but an extension of the Kane model to a crystal 
structure with lower symmetry. 



To summarize, based on the symmetry properties and 
the k • p perturbation theory, we systematically derived a 
model Hamiltonian for the 3D TI in the Bi2Ses class of 
materials. Our model Hamiltonian captures the main low 
energy physics, such as the inverted band structure and 
topologically protected surface states. The topological 
surface states have well defined spin texture, which can 
be traced back to the sign of the atomic SOC in these ma- 
terials. Furthermore, the Landau levels of a z-direction 
magnetic field for both bulk states and surface states are 
calculated. The gap of bulk Landau levels is shown to 
decrease when magnetic field increases, which may be 
observed in a magneto-optical spectroscopy. Within the 
bulk gap, the surface Landau levels appears as discrete 
peaks for the LDOS, which can be detected by STM. 
We also analyze the quantitative limitation of our model 
Hamiltonian with four bands and and show that a new 
model Hamiltonian with eight energy bands can describe 
Bi2Ses type of materials quantitatively, which will be 
useful in the future comparison with experiments. 
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Appendix A: Symmetry property of group Z)f^ 



As described in the text, the group D^^ is generated by 
a three- fold rotation operator i?3, a two- fold rotation op- 
erator R2 and an inversion operator P. It has six classes 
and correspondingly six irreducable representations, f^, 
f ^ and f^ with the upper index ± denoting the parity 
of the representation. Here we use T to denote the rep- 
resentation at r point in BZ to avoid confusion with the 
Dirac F matrices. The character table of D^^ is given in 
table (HI) [39,]. 



TABLE I: The character table for Dld{R3m) 



Dsd(3m) E 2R3 3R2 P 2PR3 3PR2 


Ft 1 1 


111 1 


f+ 1 1 


-111 -1 


f3+ 2 -1 


02-1 


fr 1 1 


1 -1 -1 -1 


f- 1 1 


-1 -1 -1 1 


f- 2 -1 


0-21 



After taking into account the spin, the C = 27r rota- 
tion induces a minus sign for the spin part, so that the 
elements of the group is doubled, which is the so-called 
double group. For D^^^ the class and irreducable repre- 
sentation of the double group is also doubled. The char- 



acter table for the double group of D^^ is given in table 
©[39]. 

When constructing the double group, it is useful to 
consider the decomposition of the direct product of ff 2 3 
and spinor represenation Fe, which is given by 



ff 






(Al) 
(A2) 
(A3) 



Furthermore when considering about the matrix ele- 
ments of k • p theory, the following direct productors will 
be quite helpful. 



(f±)*®f± = f++f+ 


+ f3+ 


(A4) 


(f+)*®f6- = fr+f2- 


+ f3 


(A5) 


(f6+r®f± = f3± 




(A6) 


{f^r®tt = tf 




(A7) 


(f+r®f± = f3± 




(A8) 


it^r^Tf = Ti 




(A9) 


(f|(5)r®f-,) = rr 




(AlO) 


(f|(5)r®f-4) = f2- 




(All) 



TABLE II: The character table for the double 
group oi DldiRSm). 



D3d{3n 


i) E 2R3 3R2 P 2PR3 3PR2 


C 2CR3 3CR2 CP 2CPR3 3CPR2 


f/ 












1 






rV 




1 -1 




1 -1 




-1 




1 -1 


rV 




-1 




-1 


2 - 


1 


2 


-1 


rV 




-1 i 




-1 i 


-1 ] 


-i 


_]^ 


1 -i 


rV 




-1 -i 




-1 -i 


-1 ; 


i 


-1 




tV 




1 




1 


-2 


1 


-2 


-1 


Ti 






-1 


-1 -1 




1 


-1 


-1 -1 


Fa 




1 -1 


-1 


-1 1 




-1 


-1 


-1 1 


Ts 




-1 


-2 


1 


2 


1 


-2 


1 


r4 




-1 i 


-1 


1 -i 


-1 1 


-i 




-1 i 


Ts 




-1 -i 


-1 


1 i 


-1 1 


i 




-1 -i 


fe" 


2 


1 


-2 


-1 


-2 - 


1 


2 


1 



Appendix B: F matrix 



l(8)r2 



l(8)r3. 



(Bl) 



The five Dirac F matrices can be defined as 

Fi = (Ji (g) Ti F2 = (72 (8) Ti F3 = (73 (g) Ti 



which satifies Clifford algebra {F^, F5} = 2Sab- The other 
ten F matrices are given by Tab = [Fa, Tij]/2i. Explicitly, 
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ab 


is 


given by 



CTi 



)Ti,(7j 



)Ti\/2i = SijkCTk ^ 1 



[(Ji (g) Ti, 1 (g) T2]/2i = CTi (g) ra 
[cr^ (g) Ti , 1 (g) r3]/2i = -cr^ g) r2 
[Ig)r2,lg)r3]/2i = lg)ri 



(B2) 
(B3) 
(B4) 
(B5) 



where z,j = 1,2,3. Now let's check the properties of 
the fifteen F matrice under the time reversal operation 
T and inversion operation P. We assume the T matrices 
are written in the basis |Plt, |), \P2~,l), |Pll,-^) 
and |P2^,— ^), then the transformation matrix of the 
symmetry operation has been obtained in Sec. IIIII With 
these transformation matrices, we have 









5^ =^-15 

In fact P operator is exactly T^ here. 



1,2,3,4 (B6) 
(B7) 



Tr46T 



-PT45P 



45- 



(B8) 

(B9) 

(BIO) 

(Bll) 



where 


^,j = 1,2,3. 






Next let's consider about R2. 








^2ri,4^2 = "^1,4, 




(B12) 




^2r2,3,5^2 = r2,3,5 




(B13) 




^2ri2, 31, 24, 34, 15, 45^2 = 


— ri2, 31, 24, 34, 15, 45 ,(214) 




^2r23,14,25,35^2 — ^23, 14, 25, 35- 


(B15) 



Finally let's talk about the three fold rotation symme- 
try. Under the rotation operation Rz{0), the T matrice 
are transformed as r'(^) = e* 2"^re~*^^, then 



dT'je) 

de 



;P,r'( 



(B16) 



Therefore, the transformation properties of F matrice un- 
der the rotation operation are determined by the commu- 
tation relation [S,F]. The commutation relations for F 
matrice are listed as follows: 



-2iFi 



(B17) 



[S,Fi] = 2iF2, [E,r2] 

[s,r3] = [E,F4] = [s,r5] = o, 

[s,ri2] = o, [s,r34]=o (bis) 

[S, Fgi] = -2ir23 p, F23] = 2ir3i (B19) 

[S, Fm] = 2iF24 P, F24] = -2iFi4 (B20) 

[S,Fi5] = 2iF25, [S,F25] = -2iFi5 (B21) 

[S,F35]=0, p,F45]=0 (B22) 

With the above commutation relations, we can easily 
solve the equantion (|B16|) and find that 

r;(l9) = Ti cos 6> -r2 sin l9. 



^2^ = Ti sin i9 + r2 cos i9 

r23W = r23COs6>-r3iSin(9, 
TsiW =r3icos6> + r23sin(9 
) = ri4 COS — r24 sin ^, 
) = ri4 sin + r24 COS 
) = ri5 COS — r25 sin ^, 
^^ ) = Fis sin + r25 cos 

^'b{^) = Ts, r34 = r34, 
r35 = r35, r45 = r45 



r'i4 
r'24(^: 
r'i5(^: 
r'25(^: 



r' 

^ 12 



= F 



12, 



(B23) 
(B24) 

(B25) 

(B26) 

(B27) 

(B28) 



The above results indicate that under the rotation R^ 
r3,4,5 and ri2,34,35,45 behave as scalars (or pseudo- 
scalars), while the three pairs of operators {r23,r3i}, 
{ri4,r24} and {ri5,r25} behave as vectors. The cor- 
responding representation for each V matrix is given in 
table (Iml). 



TABLE III: The character table of V matrice and 
the polymals of the momentum k. 



Representation T 


{ri.Fa} 


fs" 


Fs 


fr 


F4 


fj- 


Fs 


f^ + 


Fl2 


rt 


{F23,F3l} 


rt 


{Fl4,F24} 


rt 


{Fl5,F25} 


fg- + 


F34 


f+ 
^2 


F35 


fr + 


F45 


f2- + 


{k..,ky} 


fg- 


k.,kl 


fj- 


1 , K,^ ~r rvy 1 K, z 


f^ + 


\Kjx '^'1/5 ^rCxf^y j 


rt + 


l^X OKxl^y 


fr 


oKr^rvy rvy 


f2- 


\Kx + K,xK>yi l^xl^y 1 ^^J 


fg- 


{B.,By} 


f3+ 


B, 


rt 



Appendix C: Parameters in k • p theory 

In this appendix, we hope to show the detail results 
from k • p theory. First let's consider the constraint for 
the matrix elements of the momentum from the D^^ sym- 
metry. As described above, the eigen-states can be de- 
noted by |A^,a) with A = P1±,P2± and a = ±^,±|. 
The states |A, ±1/2) belong to f^ representation. For 
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|A, ±3/2), as described above, we need to re-combine 
these two states as 



Here it is more convenient of use p± = Px ^ ipy, which 
leads to 



|A±, f4) = -i=(|A±, 3/2) + |A±, -3/2)) (CI) 
|A±,f5) = ^(|A±,3/2) - |A±,-3/2)), (C2) 

which belong to V/^ and T^ representation respectively. 
The expressions ()A4p ~ ()Alip give the decomposition 
of the direct product of these states. The momentum 
Px-, Py belongs to T^ representation, while pz belongs 
to r^ representation, therefore we require that the de- 
composition of the direct product of the eigen-states 
also include T^ and F^ to obtain non-zero matrix ele- 
ments. For example, the direct product of Fg and FJ^ 

doesn't contain r2, which indicates that the matrix ele- 
ment (Ai, ±l/2|p2|A2, F4^5) is always zero. 

The symmetry operation can further help us to obtain 
the relation between different matrix elements of the mo- 
mentum. For example, due to the Rs rotation symmetry, 
we have 

(A+,-|p,|A-,--) 
= (A+,i|4i?3P.4^3|A2",-^) 
= e "" {H^2\\Px^os— -py^m—\\K^,--) 

^ (A+,i|p.|A2-,-i)=z(A+,i|p,|A2-,-i). (C3) 

Finally we can define the independent component of the 
matrix elements as follows. 

(A+,-|p,|A^,--) = (A+,--|p^|A2-,-) 

= «{Af . 2'^^'^^' "2^ " ~*^^^' "2'^'''^^' 2^ 
= ^A+.A,- (C4) 

(A+,ib,|A-,l) = -(A+,-lb,|A-,-l) 

= QA^A,- (C5) 



-i{Kt,-\py\kl,f,) = {Kt--\py\K^,f,) 



^KfAt 



(C6) 



(Af,ib,|Aj,f5)=i(Af,-ib,|Aj,f5) 



-i{A±,ib,|A^,f5) = -(A±,-ib,|A^,f5) 



^A±Af 



AT, a; 



{A+,f4bz|A-,f4)=i?A+ 

{A+,f5b,|A-,f5) = V 



(C7) 
(C8) 
(C9) 



{A+,ib+|A2-,-^) = {A+,-i|p_|A2-,^) 

= 2i^A+ A- 

{A+,ib_|A2-,-i) = (A+,-ib+|A2-,i) 

= 

{Kt,\\p.\K^,Y,) = -i{Kt,-\\p+\K^,Y,) 



2M 



A±,AT 

(Af,ib+|A2T,r4) = -i(Af,-lb_|A2^,r4) 

= 

(Af,ib_|A2^,r5) = i(Af,-lb+|A2T,r5) 



(Af,i|p+|Aj,r5) = i(A±,-ib_|Aj,r5) 
= 



(CIO) 
(Cll) 

) 

(C12) 

) 

(013) 
(C14) 

(C15) 



Q 



Time reversal symmetry indicates that P^+ ^- , 
can be chosen to be real (-P^+ ^- = P*+ .-, 



Ar.A^ 



Ar.A2 



^A+,AJ 



R 



KTAZ 



^A+,A- 
_ Q* 
■*A+,A,- 



) while M^±,T = iN*^!,T and 
. Since the matrix element between 



|Pll,±|) and |P2^,±^) is quite important, we denote 
(Plt,ib.|P2;,-i) = (Pll,-ib.|P2;,i) = 
i{Plt, \\Py\P2-+, -\) = -i{Plt, -ib.|P2;, \) 



Po 



(C16) 



(Pll,ib,|P2;,i) = -(Pll,-lb,|P2;,-i) 

= Qo (C17) 

Now we consider the perturbation theory. The degen- 
erate perturbation formulism is given by 



^mm' - ^rnA 



h: 



(1) 



m^mm' 
rr/ 



rT{2) ^\S^ TTf TTf 

-'^mm' o / V -'--'- ml -'■-'■Im' 



El 



El 



(C18) 
(C19) 



(C20) 



H. 



+ 



(3) 



^E 



l,m" 



-^^/-^/m/'-^n 



ml Im" ira"ira' 



{Era' — El) {Em" — Ei) 



^mm"-^m"l^lm' 



{Em — El) {Em" — El) 
+ "^ / .^ml^ll'^l'm' 



1 



{Em — Ei){Em — Ei>) 



1 

{Em' - El) {Em' - El') 



(C21) 



Here m and m^ are taken from |Pll,l/2) = |1), 
|P2;,l/2) = |2), |Pll,-l/2) = |3) and |P2;,-l/2) = 
1 4) with the energy Ei = Es and E2 = £^4 and Ei < E2. 
I is taken from the other bands except for these four 
bands. The expression from perturbation calculation of 
our model Hamiltonian with four bands is given as fol- 
lows and the values of the parameters are listed in table 
llVi For our new model Hamiltonian with eight band 
model, the perturbation procedure is the same to our 
model Hamiltonian with four bands and here we only list 
the values of the parameters in table |Vl 



Co + Mo = El 
Co-Mo = E2 

Ci+Mi 



2mo 



h^ 



m^ 



E 



\Q 



P1 + ,A- 



E, -E, 



A-, 1/2 



(C22) 
(C23) 

(C24) 
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2mo 



C2 + M2 = 

, |Mpi+,A-|2 

El - Ef^-^Yi 
Ci-Mi= ^' 

C2 - M2 = 



m'i 



E 



\ppi 



+ ,A- 



|A^. 



P1 + ,A- 



Ei-E 

2 



'A- 



-1/2 



\M, 



A+,P2- 



2mo 

|2 



El - ^A-,r5 

.2 Z^ 



m, 







El - £^A+,i/2 



(C25) 
(C26) 



h^ 



mf 



E 



A+,P2- 



Ei-E^ 



A+,-1/2 



I AT, 



A+,P2- 



Ei - -E'A+,r4 -El — Ef^+^r^ 

mo 

-Bo = — Qo 

Too 



For i?i and i?2 term we have 



(C27) 
(C28) 
(C29) 



ill ~ it2 = — o 



E 



|Mpi + ,A-pPpi+,P2- 

{Ep2- - Ej^-){Epi+ - Ej^ 



7 > ^Pl + Ar^A+A-^A+P2- 77^ ^^ v7^ ^^ 

Ai ,A+ \ 1 2 



) (^P2- - E^-){Ep2- - E^+) 



Ri -\- R2 = 



h^ 



m^ 



E 



Ppi+,P2- |M"p2-,A-| 

{Ep2- - £'a+)(^pi+ - ^A+) 



E Ppi^A-M^-AtK2-At ( 



1 



(C30) 



{Epi+ - E^-){Epi+ - E^+) {Ep2- - E^-){Ep2- - E^+) 
I 



(C31) 



Now we study the effect of magnetic field. Under mag- 
netic field, there are two different kinds of contribution. 
One is the orbital term, which induce the Landau level 
and has been considered in Sec. IVII The other one is 
the Zeeman type term, which is described by an effec- 
tive g factor. In the following we will discuss about 
the effective g factor in detail. There are two kinds of 
contributions to the effective g factor. One comes from 
the atomic g factor, which can be estimated from the 
ab initio calculation. In an atom, the electron spin and 
orbital angular momentum couples to magnetic field by 
Hzee = ^ {giL^gsS)-B = ^goJ-B, where J = S + L 
is the total angular momentum and go is so called Lande 
g- factor. The wave functions for the basis of our model 
Hamiltonian have been calculated from ab initio calcu- 
lation, which can be projected into the atomic orbitals. 
Since for each atomic orbitals, the g factor is simply given 

by ^0 = 1 + '^^'^^^^~2jfjX^i)^^^^^^^ ^ ^^^ effective go can 



be easily calculated, which is found to be go ^ 1.2. An- 
other contribution to the effective g factor origins from 
the second order perturbation, which is related to the cor- 
rection to the effective mass term. The relation between 
the effective mass and effective g factor in the ordinary 
semiconductors is known as the Roth's formula[4q. Here 
the second order correction to g factor is given by 



(2) 
9lz 



mo 



E 



iPpi+ 



P1 + ,A- 



|M 



P1+,A- 



Ei — E^ 

\Npi^ 



-1/2 



P1 + ,A- 



Ei 



^A-,r4 



Ei-E^ 



A-,r5 



(2) 
^1 



2) _ 4 ^ QP1+,A-^P1+,A 



£^Pi+ — E, 

|^A+,P2 



P1+ - ^A-,1/2 
2 



(C32) 
(C33) 



(2) ^ ^ V- 

^^^ mo f^ V^2- ^A+,-1/2 



TABLE IV: The summary of the parameters in 
our model Hamiltonian with four bands. 





Bi2Se3 


Bi2Te3 


Sb2Te3 


Ao(ey-A) 


3.33 


2.87 


3.40 


BoieV'A) 


2.26 


0.30 


0.84 


CoieV) 


-0.0083 


-0.18 


0.001 


C^{eV-k') 


5.74 


6.55 


-12.39 


C2{eV-A^) 


30.4 


49.68 


-10.78 


Mo{eV) 


-0.28 


-0.30 


-0.22 


Mi{eV-A^) 


6.86 


2.79 


19.64 


M2{eV'k^) 


44.5 


57.38 


48.51 


Ri{eVk^) 


50.6 


45.02 


103.20 


R2{eVA') 


-113.3 


-89.37 


-244.67 


giz 


-25.4 


-50.34 


-UA5 


dip 


-4.12 


-2.67 


-2.43 


92z 


4.10 


6.88 


14.32 


92p 


4.80 


3.43 


16.55 



|Ma+ 



P2- 



I^A+, 



P2- 



E2 



(2) 
92p 



4 



Eo 



mo 



E 

A+ 



■^A+,r5 
-Pa 



Qa+,P2-^A+,P2- 



E 



P2- 
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(C34) 
(C35) 



where giz{p) and g2z{p) are defined in (j35|). Therefore 
finally our effective g factor is the summation of the above 
two different contributions, 



ga = 9o^9^a\ 



a = Iz, 2z, Ip, 2p (C36) 



and the values of effective g factor are given in table llVi 
From table ([Y]), we find that for |Pll, ±^) band, there 
is a strong anisotropy which comes from the large con- 
tribution of the second order perturbation of the states 
|P2-,±f) and |P2I,±^). 



TABLE V: The summary of the parameters in 
the eight band effective model. 





Bi2Se3 


Bi2Te3 


Sb2Te3 


Pi{eV-A) 


3.33 


2.87 


3.40 


Qi{eV-A) 


2.26 


0.30 


0.84 


P2{eV'k) 


2.84 


2.68 


3.19 


Q2{eV-k) 


2.84 


2.68 


3.19 


P3{eV-A) 


-2.62 


-1.94 


-2.46 


Q3{eV-A) 


2.48 


1.23 


2.11 


F,{eV-A^) 


3.73 


7.16 


3.82 


Ki{eV-A^) 


6.52 


3.72 


2.49 


F3{eV-k^) 


-1.12 


3.76 


-32.03 


K3{eV-k^) 


-14.0 


-7.70 


-59.28 


F,{eV-k') 


1.50 


-0.62 


-2.26 


K,{eV-k^) 


-3.11 


-7.17 


-13.00 


Fr{eV-k') 


2.71 


3.77 


5.04 


Kr{eV-k^) 


-5.08 


22.27 


2.40 


U35 = US6{eV-k^) 


-2.31-7.45i 


-2.21-9.85i 


-ll.31-46.00i 


V35 = -V3%{eV-k^) 


-l.05-5.98i 


-2.43-3.53i 


-4.50-22.80i 


F37{eV-k^) 


2.47 


4.39 


16.96 


K37{eV-k^) 


-8.52 


-6.50 


-24.17 


U47(eV-k^) 


-7.86 


-4.29 


-45.46 


V^r{eV-k') 


-8.95i 


-0.83i 


-17.64i 


U58 = -USsieV-k^) 


-2.31-2.57i 


-0.24-3.69i 


-2.01-3.98i 


V38^VelieV-k^) 


-0.64-4.29i 


-0.85-6.64i 


1.28-9.021 


E,[eV) 


-0.29 


-0.48 


-0.22 


E3{eV) 


0.28 


0.12 


0.22 


E,{eV) 


-0.57 


-0.63 


-0.88 


Er{eV) 


-0.98 


-1.18 


-1.51 
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